緣起

土地孕育著無數生命，如同人類的母親。然而當台灣工商業愈加發達，許多農地也在發展歷程中遭汙染。

為讓民眾瞭解農地汙染狀況，台灣環境資訊協會2015年10月發起「守護農地計畫」，希望敦促政府積極行政，讓國民吃得安心，使土地生生不息。

透過智庫驅動（DSP）發起的「D4SG資料英雄計畫」，環資與政治大學資訊科學系/所、新聞系的學生2016年3月下旬組成跨領域團隊，透過資料科學方法，協助環資找出可能受重金屬汙染卻未受政府管制的農地，且分析各縣市的差異，而後希望發展相關專題報導。

這份教案在保留原汁原味的前提下，以環保署2016年的列管資料重現當時的研究成果，以本著作係採用創用 CC 姓名標示 3.0 台灣 授權條款授權。

環境設定

1.載入實作練習包

2.安裝需到的R套件

list.of.packages <- c("knitr", "dplyr", "plotly", "scales", "leaflet", "kableExtra", "stringr")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

3.資料匯入

  • 使用 read.csv 匯入資料
  • 建議把游標移到" " (雙引號) 裡面按Tab鍵選取檔案
data_EPA <- read.csv("data/環保署列管污染農地_utf8.csv", fileEncoding = "utf8")
data_TARI <- read.csv("data/農地重金屬超標未列管_utf8.csv", fileEncoding = "utf8")
# 資料瀏覽

PART1–環保署資料實作

環保署「土壤及地下水列管場址」資料，提供有公告列管的農地污染控制場址地號，截至2016-12-18共計5485筆。

1.使用 dplyr 套件製作敘述統計表

library(dplyr)
#因為沒有全部要用，擷取需要用到的資訊，並另外命名為英文
data_EPA <- data_EPA %>% 
  select(county=1, coordinate=3, area=7, control_date=8, free_date=10)

summarise_data <- data_EPA %>% group_by(county) %>% 
  summarise(count = n(),                                 # 所有案件數
            sum_area = sum(area),                        # 場址面積總和
            control_area = sum(area[free_date=="無"]),   # 列管面積總和 (平方公尺)
            free_area = sum(area[free_date!="無"]),      # 解除面積總和 (平方公尺)
            count_control = sum(free_date=="無"),        # 列管案件
            count_free = sum(free_date!="無")) %>%       # 解除列管案件
  mutate(avg_area = sum_area/count,                      # 平均場址面積 (平方公尺)
         ratio_control = sprintf("%1.0f%%",count_control/count*100), # 列管比例
         ratio_free = sprintf("%1.0f%%",count_free/count*100))       # 解除列管比例

2.計算平均列管月份

data_EPA <- mutate(data_EPA,
                   control_date = as.Date(control_date), # 把列管時間轉成日期格式
                   free_date = as.Date(free_date), # 把列管時間轉成日期格式
                   # 把尚未解除列管的時間帶資料收集截止時間"2016-12-18"，以便計算列管時間
                   free_date = replace(free_date, 
                                       which(is.na(free_date)), 
                                       as.Date("2016-12-18"))
                   )

# 計算列管月份
data_EPA$totol_month <- NA
for(i in 1:nrow(data_EPA)){
  tryCatch({ # 因為資料裡有其中一筆沒有列管時間，我們加tryCatch讓他忽略錯誤，讓迴圈可以正常執行完
  data_EPA$totol_month[i] <- 
    length(seq(from=data_EPA$control_date[i], to=data_EPA$free_date[i], by='month')) 
  }, error=function(e){})
}

# 計算各縣市列管中與解除列管的平均列管月份
avg_month <- data_EPA %>% group_by(county) %>% 
  summarise(control_month = mean(totol_month[free_date=="2016-12-18"], na.rm = TRUE) %>% 
                            round(), 
            free_month = mean(totol_month[free_date!="2016-12-18"], na.rm = TRUE) %>% 
                         round())
# 因為資料中有一筆沒有列管時間，因此把參數 na.rm 改為TRUE，才可以運算

# 把無列管中的平均月份帶0
avg_month[is.na(avg_month)] <- 0

# 合併回summarise_data
summarise_data <- summarise_data %>% 
  left_join(avg_month, by="county")

3.計算全國加總

tmp <- colSums(summarise_data[,2:8]) %>% t() %>% as.data.frame() %>% 
  mutate(ratio_control = sprintf("%1.0f%%",count_control/count*100),
         ratio_free = sprintf("%1.0f%%",count_free/count*100)) %>% 
  cbind(.,colMeans(summarise_data[,11:12]) %>% round() %>% t() %>% as.data.frame()) %>%
  mutate(county="全國") %>% 
  select(12,1:11)

summarise_data <- rbind(tmp, summarise_data)

4.將欄位改成中文敘述並進行展示

tmp2 <- 
  summarise_data %>% 
  select(`全台縣市`=county,
         `案件數`=count,
         `場址面積總和`=sum_area,
         `列管面積總和`=control_area,
         `解除列管面積總和`=free_area,
         `平均場址面積`=avg_area,
         `列管案件`=count_control,
         `解除列管案件`=count_free,
         `列管案件比`=ratio_control,
         `解除列管案件比`=ratio_free,
         `列管中_平均列管月份`=control_month,
         `解除列管_平均列管月份`=free_month) %>% 
  arrange(desc(`案件數`))
tmp2
全台縣市 案件數 場址面積總和 列管面積總和 解除列管面積總和 平均場址面積 列管案件 解除列管案件 列管案件比 解除列管案件比 列管中_平均列管月份 解除列管_平均列管月份
全國 5485 9263582.55 3576260.23 5687322.33 51145.5310 2468 3017 45% 55% 40 35
彰化縣 2503 5038536.34 2282395.32 2756141.02 2012.9989 1340 1163 54% 46% 25 28
桃園市 1890 2511211.33 1216032.30 1295179.03 1328.6832 1060 830 56% 44% 48 42
臺中市 600 742225.90 10292.02 731933.88 1237.0432 11 589 2% 98% 32 30
新竹市 202 362010.73 2746.50 359264.23 1792.1323 2 200 1% 99% 7 27
臺南市 104 194793.93 34659.06 160134.87 1873.0186 37 67 36% 64% 105 30
高雄市 49 84948.44 0.00 84948.44 1733.6416 0 49 0% 100% 0 32
苗栗縣 37 43766.63 5729.23 38037.40 1182.8819 5 32 14% 86% 12 27
雲林縣 25 58230.54 4329.00 53901.54 2329.2216 2 23 8% 92% 78 31
臺北市 22 48852.15 0.00 48852.15 2220.5523 0 22 0% 100% 0 37
嘉義市 19 46035.64 12938.00 33097.64 2422.9284 5 14 26% 74% 28 49
新北市 13 37235.00 0.00 37235.00 2864.2308 0 13 0% 100% 0 22
南投縣 11 4046.00 2747.00 1299.00 367.8182 4 7 36% 64% 167 103
宜蘭縣 5 11798.15 3156.64 8641.51 2359.6300 1 4 20% 80% 98 22
屏東縣 3 75150.81 1235.16 73915.65 25050.2700 1 2 33% 67% 6 22
嘉義縣 2 4740.96 0.00 4740.96 2370.4800 0 2 0% 100% 0 16

5.使用 plotly 套件繪製統計圖表

  • Bar Chart
library(plotly)

plot_ly(summarise_data, x = ~county, y = ~count_control, type = "bar") %>% 
  layout(title = "所有案件",
         xaxis = list(title = '全台縣市'),
         yaxis = list(title = '案件數'),
         width = 750, height = 400)
# x軸依據案件數做排序
plot_ly(summarise_data, 
         x = ~reorder(county, -count_control), 
         y = ~count_control, type = "bar") %>% 
  layout(title = "所有案件",
         xaxis = list(title = '全台縣市'),
         yaxis = list(title = '案件數'),
         width = 750, height = 400)
  • Grouped Bar Chart
  • with Hover Text and Rotated Labels
plot_ly(summarise_data, 
        x = ~reorder(county, -count_control), 
        y = ~count_control, type = 'bar',
        name = '列管案件', text = ~ratio_control) %>%
  add_trace(y = ~count_free, name = '解除列管案件', text = ~ratio_free) %>%
  layout(xaxis = list(title = "", tickangle = -45), 
         yaxis = list(title = 'count'), barmode = 'group',
         width = 800, height = 400)
  • Pie Chart
#列管縣市太多，很多面積很小，取場址面積大於中位數的來畫圖
summarise_data1 <- summarise_data %>% 
  arrange(sum_area %>% desc) %>%  # 降冪排序
  filter(sum_area > median(sum_area)) %>%  # 篩選大於中位數者
  slice(-1) # 移除全國總和數據

plot_ly(summarise_data1, labels = ~county, values = ~sum_area, type = 'pie',
        textinfo = 'label+percent') %>%
  layout(title = '場址面積總和',  width = 800, height = 500,
         margin=list(l = 100, r = 50, b = 100, t = 100, pad = 4))

6.場址座標轉經緯度

  • 把TWD97位址資訊抓出來，分割成x軸與y軸
head(data_EPA,3) #可以看到座標欄位裡面有：和，
  county            coordinate area control_date  free_date totol_month
1 臺北市 X：301277，Y：2777431  387   2002-10-04 2004-11-02          25
2 臺北市 X：301267，Y：2777372 7738   2002-10-04 2004-11-02          25
3 臺北市 X：301310，Y：2777998  580   2004-01-27 2005-03-01          14
TWD97 <- data_EPA$coordinate %>% as.character() %>% 
  stringr::str_split(.,'[，：]',simplify = TRUE) %>% 
  as.data.frame(., stringsAsFactors=FALSE) %>% select(x=2,y=4) 
source("https://raw.githubusercontent.com/snexuz/TWD97TM2toWGS84/master/TWD97TM2toWGS84.R")
WGS84 <- TWD97TM2toWGS84(TWD97$x, TWD97$y) %>% as.data.frame()

#合併回原本的data
data_EPA$coordinate <- NULL
data_EPA <- cbind(data_EPA,TWD97,WGS84)
head(data_EPA)
  county    area control_date  free_date totol_month      x       y      lat      lon
1 臺北市  387.00   2002-10-04 2004-11-02          25 301277 2777431 25.10434 121.5084
2 臺北市 7738.00   2002-10-04 2004-11-02          25 301267 2777372 25.10381 121.5083
3 臺北市  580.00   2004-01-27 2005-03-01          14 301310 2777998 25.10946 121.5088
4 臺北市 3205.65   2004-12-10 2008-04-18          41 299512 2779782 25.12562 121.4910
5 臺北市 2855.14   2004-12-10 2008-04-18          41 299518 2779865 25.12637 121.4911
6 臺北市 3513.00   2004-12-10 2008-04-18          41 299669 2779961 25.12723 121.4926

PART2–農試所資料實作

農委會農業試驗所（農試所）之「土壤品質及生產力調查」資料。原始資料是1992年至2008年間進行的全國土壤採樣調查資料，累積約13萬筆，涵蓋78萬公頃表土，包含鎘、鉻、銅、鎳、鉛、鋅等六項重金屬有效性之調查。該資料是以250公尺*250公尺網格（6.25公頃）為單位。

由於農試所資料與環保署現今重金屬管制標準檢測方法不同，前者是透過0.1M鹽酸方法，後者是經過王水消化法，需要進行公式轉換，在此提供933筆環保數超標數據。

1.重金屬轉換公式

農試所以 0.1M 鹽酸萃取抽出分析土壤中的重金屬濃度，其數據轉為王水消化法的推估值 (AR，以下為農試所使用的迴歸轉換公式) 後，再利用環保署現行管制標準進行篩檢。

重金屬名稱 重金屬轉換公式 環保署農地管制值
銅 CUAR = 2.035*CU + 11.884 200
鋅 ZNAR = 2.487*ZN + 89.711 600
鎘 CDAR = 1.4578*CD + 0.0323 5
鎘 CRAR = 17.35*CR+ 31.91 250
鎳 NIAR = 5.13*NI+ 14.56 200
鉛 PBAR = 2.811*PB+ 6.715 500 
data_TARI <- data_TARI %>% 
  mutate(CUAR = 2.035*CU + 11.884,
         CUAR_OVER = ifelse(CUAR>200,1,0),
         ZNAR = 2.487*ZN + 89.711,
         ZNAR_OVER = ifelse(ZNAR>600,1,0),
         CDAR = 1.4578*CD + 0.0323,
         CDAR_OVER = ifelse(CDAR>5,1,0),
         CRAR = 17.35*CR+ 31.91,
         CRAR_OVER = ifelse(CRAR>250,1,0),
         NIAR = 5.13*NI+ 14.56,
         NIAR_OVER = ifelse(NIAR>200,1,0),
         PBAR = 2.811*PB+ 6.715,
         PBAR_OVER = ifelse(PBAR>500,1,0),
         OVER = ifelse(CUAR_OVER+ZNAR_OVER+CDAR_OVER+CRAR_OVER+
                       NIAR_OVER+PBAR_OVER>=1, 1, 0))

2.計算農地樣區是否有超標情形

  • mindistance (與最近的環保署列管農地距離，單位公尺)
  • Is_monitored (境內是否有環保署列管農地。1有，0無)
  • minover (方圓1公里內其他超標點總數)
# 篩選出截至 2016-12-18 環保署仍未解除列管的農地
data_EPA_monitor <- data_EPA %>% 
  filter(free_date=="2016-12-18")

# 以農試所樣區資料為主體，計算每一個樣區與環保署列管農地的最短距離
for(i in 1:nrow(data_TARI)){
  data_TARI$mindistance[i] <- 
    ((data_TARI$XLO[i]-as.integer(data_EPA_monitor$x))^2+
       (data_TARI$YLO[i]-as.integer(data_EPA_monitor$y))^2) %>% 
    sqrt %>% min(na.rm = TRUE)
}

# 以農試所樣區為主體，計算每一個樣區方圓一公里內有多少個超標樣區
for(i in 1:nrow(data_TARI)){
  data_TARI$minover[i] <-
    ((data_TARI$XLO[i]-data_TARI$XLO)^2+(data_TARI$YLO[i]-data_TARI$YLO)^2) %>% 
    sqrt %>% sum(na.rm = TRUE)
}

# 以農試所樣區為主體，計算每一個樣區境內是否有環保署列管農地
data_TARI <- data_TARI %>% mutate(Is_monitored=ifelse(mindistance<=200, 1, 0))

3.與環保署的資料合併，計算統計指標

library(scales)
stat_over <- data_TARI %>% group_by(county) %>%
  summarise(`重金屬超標`=sum(OVER),
            `超標未列管`=sum(ifelse(OVER+Is_monitored==1,1,0))) %>%
  mutate(`超標未列管比例`=percent(`超標未列管`/`重金屬超標`),
         `超標列管比例`=percent(1-`超標未列管`/`重金屬超標`)) %>% 
  arrange(desc(`重金屬超標`))
# stat_over
county 重金屬超標 超標未列管 超標未列管比例 超標列管比例
桃園市 411 379 92.2% 7.8%
彰化縣 159 107 67.3% 32.7%
屏東縣 86 86 100.0% 0.0%
高雄市 85 85 100.0% 0.0%
台中市 77 77 100.0% 0.0%
台南市 32 32 100.0% 0.0%
新竹縣 30 30 100.0% 0.0%
嘉義縣 23 23 100.0% 0.0%
宜蘭縣 22 22 100.0% 0.0%
苗栗縣 18 17 94.4% 5.6%
南投縣 13 13 100.0% 0.0%
新北市 11 11 100.0% 0.0%
雲林縣 10 10 100.0% 0.0%
花蓮縣 8 8 100.0% 0.0%
台東縣 5 5 100.0% 0.0%
嘉義市 2 2 100.0% 0.0%
台北市 1 1 100.0% 0.0%

PART3–用leaflet畫地圖

地圖一：兩份資料交叉比對

library(leaflet)

#先把要點上圖的資料篩選出來
EPA <- data_EPA %>% filter(free_date=="2016-12-18") %>% select(lat, lon, area)
TARI <- data_TARI %>% select(lat = Y,lon = X, area = AREA)


leaflet() %>% setView(lng=120.58,lat=23.58,zoom=8) %>% 
  addProviderTiles("Esri.WorldImagery") %>%
  addCircles(data = EPA, color = "red", 
             lng = ~lon, lat = ~lat, weight = 1, radius = ~sqrt(area/2)) %>%
  addCircles(data = TARI, color = "yellow",
             lng = ~lon, lat = ~lat, weight = 1, radius = ~sqrt(area/2))

地圖二：用log(mindistance)畫地圖

TARI <- data_TARI %>% 
  select(lat = Y,lon = X, area = AREA ,mindistance) %>%
  mutate(log_mindistance=log(mindistance))

leaflet() %>% 
  addCircles(data=TARI,lng = ~lon, lat = ~lat, 
             radius = ~sqrt(area/2), weight = 1,
             fill=TRUE, fillOpacity = 0.8,
             color=~colorNumeric(c("#CD0000", "#FFFFFF","#0B752F"),
                                 log_mindistance)(log_mindistance)) %>% 
  addProviderTiles(providers$CartoDB.Positron) %>% 
  addLegend(position = 'topleft',
            pal =colorNumeric(c("#CD0000", "#FFFFFF","#0B752F"),
                              domain=TARI$log_mindistance),
            values=TARI$log_mindistance,
            title = 'log-mindistance')